! -----------------------------------------------------------------------
! File:             matrix_operators.f90
! Purpose:          Definition of vectorized operators and related operations
! Contains:
! Revision History: Ver. 1.0 May. 2007 Alberto Scotti
! -----------------------------------------------------------------------

!****************
MODULE Matrix_Oper
!****************
!Definitions of data structures associated to the problem
USE ElliDef 
! Metric Coefficients
USE Metric, ONLY: g,g_mask
! basic matrix operators
USE Basic_matrix

IMPLICIT NONE
SAVE 
PRIVATE

! Public interfaces
PUBLIC GRAD,adv,operation,FREEUPCRS,FREEUPDIA,FREEUP,ADB, &
     DIAGDA,gijdj,EXTRACT_U_A,EXTRACT_U_P,FILTER_Z_A,FILTER_Z_P,adv_sym, &
     EXTRAPOLATE_BC,APB,CRS2DIA

INTEGER :: ierr
INTEGER, PARAMETER :: errOK=0

!REAL (KIND=Precision), PARAMETER :: zero=0.0D0
! GLOBALLY AVAILABLE DEFITIONS (WITHIN THIS MODULE)
INTERFACE
   subroutine diamua (nrow,job, a, ja, ia, diag, b, jb, ib)
      real*8 a(*), b(*), diag(nrow) 
      integer ja(*),jb(*), ia(nrow+1),ib(nrow+1) ,nrow,job
    END subroutine diamua
 END INTERFACE

!!$INTERFACE EXTRACT_U
!!$   MODULE PROCEDURE EXTRACT_U_A,EXTRACT_U_P
!!$END INTERFACE
!!$
!!$INTERFACE FILTER_Z
!!$   MODULE PROCEDURE FILTER_Z_A,FILTER_Z_P
!!$END INTERFACE
INTERFACE FREEUP
   MODULE PROCEDURE FREEUPCRS,FREEUPDIA,FREEUPM_P
END INTERFACE

INTERFACE APB
   MODULE PROCEDURE APB_CRS,APB_DIA,APB_MP
END INTERFACE

INTERFACE REALLOCATE
   MODULE PROCEDURE REALLOCATE_CRS
END INTERFACE


CONTAINS


 FUNCTION GRAD(Info,i)
! Returns pointers to Allocated and filled  
! compact form of the vectorized gradient operator along dir i
! based on DR boundary conditions along all sides. It applies to a scalar field
! defined on (0:nx+1,0:ny+1,0:nz+1). On the ghost points DR conditions are 
! expected. IF NM conditions are given, DR conditions have to be extrapolated.
! interface declaration
TYPE (Elli) :: Info
TYPE (CRS), POINTER :: GRAD
TYPE (CRS), POINTER :: A=>NULL(),B=>NULL(),C=>NULL()

INTEGER :: i,n





   A=>gijdj(Info,i,1,'b')

   B=>gijdj(Info,i,2,'b')
   IF(Info%nDim==2) THEN
      GRAD=>APB(A,B)
      CALL FREEUPCRS(A);CALL FREEUPCRS(B)
      RETURN
   ELSE
      C=>APB(A,B) ! C=A+B

      CALL FREEUPCRS(A);CALL FREEUPCRS(B)

      A=>gijdj(Info,i,3,'b')

      GRAD=>APB(A,C) ! Info%G(i)%crsptr=A+C

      CALL FREEUPCRS(A);CALL FREEUPCRS(C)


   ENDIF

END FUNCTION GRAD


   


FUNCTION EXTRACT_U_A(G,keep) RESULT (EXTRACT_U)
! extract the upper triangular part from a CRS represented matrix
! NOTE IT IS DONE IN PLACE, IT IS MEANT TO BE INVOKED AS
! G=>EXTRACT_U(G)
TYPE(CRS), POINTER :: EXTRACT_U
TYPE(CRS) :: G
INTEGER :: n
LOGICAL, OPTIONAL :: keep
ALLOCATE(EXTRACT_U)

n=SIZE(G%rowIndex)-1 ! number of rows

CALL getu(n,G%rValues,G%columns,G%rowIndex,G%rValues,G%columns,G%rowIndex)

ALLOCATE(EXTRACT_U%rValues(G%rowIndex(n+1)-1),EXTRACT_U%columns(G%rowIndex(n+1)-1),EXTRACT_U%rowIndex(n+1))

EXTRACT_U%rValues=G%rValues(1:G%rowIndex(n+1)-1)
EXTRACT_U%columns=G%columns(1:G%rowIndex(n+1)-1)
EXTRACT_U%rowIndex=G%rowIndex
EXTRACT_U%nNonZeros=G%nNonZeros
IF((.NOT.PRESENT(keep)).OR.(.NOT.keep)) &
     DEALLOCATE(G%rValues,G%columns,G%rowIndex)



END FUNCTION EXTRACT_U_A


FUNCTION EXTRACT_U_P(G,keep) RESULT (EXTRACT_U)
! extract the upper triangular part from a CRS represented matrix
! NOTE IT IS DONE IN PLACE, IT IS MEANT TO BE INVOKED AS
! G=>EXTRACT_U(G)
TYPE(CRS), POINTER :: EXTRACT_U
TYPE(CRS), POINTER :: G
INTEGER :: n
LOGICAL, OPTIONAL :: keep
ALLOCATE(EXTRACT_U)

n=SIZE(G%rowIndex)-1 ! number of rows

CALL getu(n,G%rValues,G%columns,G%rowIndex,G%rValues,G%columns,G%rowIndex)

ALLOCATE(EXTRACT_U%rValues(G%rowIndex(n+1)-1),EXTRACT_U%columns(G%rowIndex(n+1)-1),EXTRACT_U%rowIndex(n+1))

EXTRACT_U%rValues=G%rValues(1:G%rowIndex(n+1)-1)
EXTRACT_U%columns=G%columns(1:G%rowIndex(n+1)-1)
EXTRACT_U%rowIndex=G%rowIndex
EXTRACT_U%nNonZeros=G%nNonZeros
IF((.NOT.PRESENT(keep)).OR.(.NOT.keep)) &
     CALL FREEUPCRS(G)



END FUNCTION EXTRACT_U_P


FUNCTION FILTER_Z_A(G,tol,keep) RESULT (FILTER_Z)
! removes entries less than tol from a CRS represented matrix
! NOTE IT IS DONE IN PLACE, IT IS MEANT TO BE INVOKED AS
! G=>EXTRACT_U(G)
INTERFACE
   SUBROUTINE filter(n,job,tol,a,ja,ia,b,jb,ib,len,ierr)
     integer n,job,ja(*),ia(*),jb(*),ib(*),len,ierr
     real*8 a(*),b(*),tol
   END SUBROUTINE filter
END INTERFACE
TYPE(CRS), POINTER :: FILTER_Z
TYPE(CRS) :: G
REAL (KIND=Precision) :: tol
INTEGER :: n,ierr
LOGICAL, OPTIONAL :: keep
ALLOCATE(FILTER_Z)
ierr=0 ! bug with filter, ierr is not defined bu filter
! below if everything is all right...
n=SIZE(G%rowIndex)-1 ! number of rows

CALL filter(n,1,tol,G%rValues,G%columns,G%rowIndex, &
     G%rValues,G%columns,G%rowIndex, &
     SIZE(G%rValues),ierr)
IF(ierr/=0) THEN
   PRINT*,'filter operation has encountered an error ',ierr
   STOP
ENDIF

ALLOCATE(FILTER_Z%rValues(G%rowIndex(n+1)-1),FILTER_Z%columns(G%rowIndex(n+1)-1),FILTER_Z%rowIndex(n+1))
FILTER_Z%rValues=G%rValues(1:G%rowIndex(n+1)-1)
FILTER_Z%columns=G%columns(1:G%rowIndex(n+1)-1)
FILTER_Z%rowIndex=G%rowIndex
FILTER_Z%nNonZeros=G%nNonzeros
IF((.NOT.PRESENT(keep)).OR.(.NOT.keep)) &
     DEALLOCATE(G%rValues,G%columns,G%rowIndex)
! CAN'T USE FREEUPCRS, G is not a pointer


END FUNCTION FILTER_Z_A


FUNCTION FILTER_Z_P(G,tol,keep) RESULT (FILTER_Z)
! removes entries less than tol from a CRS represented matrix
! NOTE IT IS DONE IN PLACE, IT IS MEANT TO BE INVOKED AS
! G=>EXTRACT_U(G)
INTERFACE
   SUBROUTINE filter(n,job,tol,a,ja,ia,b,jb,ib,len,ierr)
     integer n,job,ja(*),ia(*),jb(*),ib(*),len,ierr
     real*8 a(*),b(*),tol
   END SUBROUTINE filter
END INTERFACE
TYPE(CRS), POINTER :: FILTER_Z
TYPE(CRS), POINTER :: G
REAL (KIND=Precision) :: tol
INTEGER :: n,ierr
LOGICAL, OPTIONAL :: keep
ALLOCATE(FILTER_Z)
ierr=0 ! bug with filter, ierr is not defined bu filter
! below if everything is all right...
n=SIZE(G%rowIndex)-1 ! number of rows

CALL filter(n,1,tol,G%rValues,G%columns,G%rowIndex, &
     G%rValues,G%columns,G%rowIndex, &
     SIZE(G%rValues),ierr)
IF(ierr/=0) THEN
   PRINT*,'filter operation has encountered an error ',ierr
   STOP
ENDIF

ALLOCATE(FILTER_Z%rValues(G%rowIndex(n+1)-1),FILTER_Z%columns(G%rowIndex(n+1)-1),FILTER_Z%rowIndex(n+1))
FILTER_Z%rValues=G%rValues(1:G%rowIndex(n+1)-1)
FILTER_Z%columns=G%columns(1:G%rowIndex(n+1)-1)
FILTER_Z%rowIndex=G%rowIndex
FILTER_Z%nNonZeros=G%nNonzeros
IF((.NOT.PRESENT(keep)).OR.(.NOT.keep)) &
     CALL FREEUPCRS(G)

END FUNCTION FILTER_Z_P



 
FUNCTION adv(G,v,di)

! object oriented sparse matrix-vector product.
! hides the details of the operation from the user
! matrix represented in CSR format
REAL (KIND=Precision), POINTER, DIMENSION (:,:,:) :: adv
TYPE (CRS)  :: G
REAL (KIND=Precision), INTENT(IN),  DIMENSION (:,:,:) :: v
INTEGER, INTENT (IN), DIMENSION (3) :: di
INTEGER :: nx,ny,nz


nx=SIZE(v,1)+di(1);ny=SIZE(v,2)+di(2);nz=SIZE(v,3)+di(3)

   ALLOCATE(adv(nx,ny,nz))


CALL mkl_dcsrgemv('n',SIZE(G%rowIndex)-1,G%rValues,G%rowIndex,G%columns,v,adv)


END FUNCTION adv

FUNCTION adv_sym(G,v,di)

! object oriented sparse matrix-vector product.
! hides the details of the operation from the user
! matrix represented in CSR format
REAL (KIND=Precision), POINTER, DIMENSION (:,:,:) :: adv_sym
TYPE (CRS)  :: G
REAL (KIND=Precision), INTENT(IN),  DIMENSION (:,:,:) :: v
INTEGER, INTENT (IN), DIMENSION (3) :: di
INTEGER :: nx,ny,nz

!!$IF(.NOT.(ASSOCIATED(v))) THEN
!!$   PRINT*,'Vector not allocated'
!!$   PRINT*,'exiting now'
!!$   STOP
!!$ENDIF
!!$IF(ASSOCIATED(adv)) THEN
!!$   DEALLOCATE(adv);NULLIFY(adv)
!!$ENDIF
nx=SIZE(v,1)+di(1);ny=SIZE(v,2)+di(2);nz=SIZE(v,3)+di(3)

   ALLOCATE(adv_sym(nx,ny,nz))


CALL mkl_dcsrsymv('U',SIZE(G%rowIndex)-1,G%rValues,G%rowIndex,G%columns,v,adv_sym)


END FUNCTION adv_sym

FUNCTION operation(type,p,dir,delta)
! creates the matrix corrisponding to operator type acting along direction dir
! on matrix p. Note that p is only used to infer the dimension of the matrix
! it does not need to be initialised after allocation.
! ('a' averaging, 'd' is difference). 
! For differencing operators delta is the grid spacing
TYPE (CRS),POINTER :: operation

INTEGER, INTENT(IN) :: dir
REAL (KIND=Precision), POINTER, DIMENSION (:,:,:) :: p
REAL (KIND=Precision), OPTIONAL, INTENT (IN) :: delta
CHARACTER(LEN=5), INTENT (IN) :: type

INTEGER :: nx,ny,nz


IF(dir<0.OR.dir>3) THEN
   PRINT*,'error operation'
   PRINT*,'dir is out of range'
   STOP
ENDIF
ALLOCATE(operation) ! this is required to allocate the entries in the derived 
                    ! data type
nx=SIZE(p,1);ny=SIZE(p,2);nz=SIZE(p,3)
SELECT CASE (TRIM(type))

CASE ('gc2eb')

   IF(.NOT.PRESENT(delta)) THEN
      PRINT*,'gc2eb operations need delta'
      STOP
   ENDIF
   
   CALL gc2eb(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,delta,dir,nx,ny,nz)
   operation%nNonzeros=nx*ny*nz ! # of columns

CASE ('c2eb')
   
   CALL c2eb(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,dir,nx,ny,nz)
   operation%nNonzeros=nx*ny*nz ! # of columns
CASE ('e2cs')
   IF(dir==1) THEN ! hack given how e2c is defined
      nx=nx-1
   ELSEIF(dir==2) THEN
      ny=ny-1
   ELSEIF(dir==3) THEN
      nz=nz-1
   ENDIF
   
   CALL e2c(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,dir,nx,ny,nz)
   
   IF(dir==1) THEN ! hack given how e2c is defined
      nx=nx+1
   ELSEIF(dir==2) THEN
      ny=ny+1
   ELSEIF(dir==3) THEN
      nz=nz+1
   ENDIF
   operation%nNonZeros=nx*ny*nz
   
CASE ('ge2cs')
   
   IF(.NOT.PRESENT(delta)) THEN
      PRINT*,'gc2es operations need delta'
      STOP
   ENDIF
   
   IF(dir==1) THEN ! hack given how ge2c is defined
      nx=nx-1
   ELSEIF(dir==2) THEN
      ny=ny-1
   ELSEIF(dir==3) THEN
      nz=nz-1
   ENDIF
   
   CALL ge2c(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,delta,dir,nx,ny,nz)
   
   IF(dir==1) THEN ! hack given how e2c is defined
      nx=nx+1
   ELSEIF(dir==2) THEN
      ny=ny+1
   ELSEIF(dir==3) THEN
      nz=nz+1
   ENDIF
   operation%nNonZeros=nx*ny*nz
   
CASE ('strip')
! remove ghost layers along direction dir
      CALL strip(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,dir,nx,ny,nz)
   operation%nNonZeros=nx*ny*nz


CASE ('zeros')
      CALL zeroes(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,dir,nx,ny,nz)
   operation%nNonZeros=nx*ny*nz

CASE ('zeroi')
      CALL zeroi(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,dir,nx,ny,nz)
   operation%nNonZeros=nx*ny*nz

CASE ('dress')
! add ghost layers along dir. If dir=0 ghost layers are added along all three directions
   CALL dress(operation%rValues,operation%columns,operation%rowIndex,operation%nNonzeros,dir,nx,ny,nz)
   operation%nNonZeros=nx*ny*nz

CASE DEFAULT
   PRINT*,type,' OPERATION NOT IMPLEMENTED YET'
   PRINT*,'in operation'
   STOP


END SELECT

END FUNCTION operation

FUNCTION ADB(C,B)
  
  ! returns in ADB the CSR representation of the product 
  ! of two CSR sparse matrixes C and B ADB=C*B 
  TYPE (CRS), POINTER :: ADB
  TYPE (CRS),POINTER :: C,B


  ! local storage
  REAL (KIND=Precision), ALLOCATABLE, DIMENSION (:) :: a
  INTEGER, ALLOCATABLE, DIMENSION (:) :: ia,ja,iw
  INTEGER :: ncol,nrow,ndummy,ierr

  IF((C%nNonZeros/=SIZE(B%rowIndex)-1)) THEN
     PRINT*,'Matrixes must identical inner dimensions'
     PRINT*,'to be multiplied'
     PRINT*,'row of first',C%nNonZeros
     PRINT*,'column of second',SIZE(B%rowIndex)-1
     
     PRINT*,'Fix and restart'
     STOP
  ENDIF


  ndummy=B%nNonZeros*20
  nrow=size(C%rowIndex)-1 ! rows in A

  ncol=B%nNonZeros
  ierr=1

  ALLOCATE(iw(ncol),ia(size(C%rowIndex)))
  DO WHILE(ierr/=0)

     ALLOCATE(A(ndummy),ja(ndummy))

     CALL AMUB(nrow,ncol,1,C%rValues,C%columns,C%rowIndex,B%rValues,B%columns,B%rowINdex,a,ja,ia,ndummy,iw,ierr)

     IF(ierr/=0) THEN
        DEALLOCATE(A,ja)
        ndummy=ndummy +B%nNonZeros*10! increase the size until right...
     
     ENDIF

  ENDDO
  
  ALLOCATE(ADB)
  ADB%nNonZeros=B%nNonZeros
!  dummy%nNonZeros=B%nNonZeros
  ndummy=IA(size(IA))-1;
!!$
  ALLOCATE(ADB%rValues(ndummy),ADB%columns(ndummy))
  ALLOCATE(ADB%rowIndex(SIZE(IA)))

  ADB%rowIndex=IA
  ADB%rvalues=A(1:IA(SIZE(IA))-1);
  ADB%columns=JA(1:IA(SIZE(IA))-1);

  
!!$
  

  DEALLOCATE(iw,A,ja,ia)

END FUNCTION ADB
FUNCTION APB_DIA(A,B,freeA,freeB) RESULT(APB)
! Return a pointer (APB) to the DIA representation of the
! sum of two DIA sparse matrices A and B APB=A+B
TYPE (DIA), POINTER :: APB
TYPE (DIA), POINTER :: A,B
LOGICAL, INTENT (IN), OPTIONAL :: freeA,freeB
! local storage
INTEGER, DIMENSION (:), ALLOCATABLE :: dist1
INTEGER, DIMENSION (:), ALLOCATABLE :: dist2
INTEGER, DIMENSION (:), ALLOCATABLE :: distance
INTEGER, DIMENSION (:,:), ALLOCATABLE::dummy

INTEGER :: p,m1,m2,m,mmax,i,ncol,nrows

! 1.0 determine if one or more input is Zero (null pointer)
IF(.NOT.ASSOCIATED(A)) THEN! A is zero
   IF(.NOT.ASSOCIATED(B)) THEN
      APB=>NULL() ! sum is zero
      RETURN
   ELSE
      ALLOCATE(APB)
      ALLOCATE(APB%rValues(SIZE(B%rValues,1),SIZE(B%rValues,2)))
      APB%rValues=B%rValues
      ALLOCATE(APB%distance(SIZE(B%distance)))
      APB%distance=B%distance
      APB%ncol=B%ncol

      IF(PRESENT(freeB).AND.freeB) CALL FREEUP(B)
      RETURN
   END IF
END IF
!
IF(.NOT.ASSOCIATED(B)) THEN! A is zero
   IF(.NOT.ASSOCIATED(A)) THEN
      APB=>NULL() ! sum is zero
      RETURN
   ELSE
      ALLOCATE(APB)
      ALLOCATE(APB%rValues(SIZE(A%rValues,1),SIZE(A%rValues,2)))
      APB%rValues=A%rValues
      ALLOCATE(APB%distance(size(A%distance)))
      APB%distance=A%distance
      APB%ncol=A%ncol

      IF(PRESENT(freeA).AND.freeA) CALL FREEUP(A)
      RETURN
   END IF
END IF
! 0.0 sanity check

IF((SIZE(A%rValues,1)/=SIZE(B%rValues,1)).OR.(A%ncol/=B%ncol)) THEN
   PRINT*,'APB_DIA cannot add matrices of different sizes'
   STOP
END IF

! 1.1 both matrixes are not zero. get to work
ALLOCATE(dist1(SIZE(A%distance)),dist2(SIZE(B%distance)))
ALLOCATE(distance(SIZE(dist1)+SIZE(dist2)))
ALLOCATE(dummy(SIZE(distance),0:2))
! First, determine the structure of the sum matrix
dist1=A%distance;dist2=B%distance
ncol=B%ncol;nrows=SIZE(A%rValues,1)
p=0;mmax=MAX(MAXVAL(dist1),MAXVAL(dist2))
DO ! parse the diagonals of A and B and determine which are present
   ! and which overlap
   m1=MINVAL(dist1);m2=MINVAL(dist2)
   m=MIN(m1,m2);IF(m>mmax) EXIT
   p=p+1;distance(p)=m
   SELECT CASE(m1-m2)
   CASE(0) ! both matrices have the diagonal
      dummy(p,0)=0 ! flag both columns need to be added

      dummy(p,1:1)=MINLOC(dist1) ! store column which holds diagonal

      dummy(p,2:2)=MINLOC(dist2)
      dist1(MINLOC(dist1))=mmax+1 ! move it away
      dist2(MINLOC(dist2))=mmax+1

   CASE (1:) ! diagonal comes from B
      dummy(p,0)=2
      
      dummy(p,2:2)=MINLOC(dist2)
      dist2(MINLOC(dist2))=mmax+1

   CASE(:-1) ! diagonal comes from A
      dummy(p,0)=1
      
      dummy(p,1:1)=MINLOC(dist1)
      dist1(MINLOC(dist1))=mmax+1

   END SELECT

END DO
ALLOCATE(APB)
APB%ncol=ncol;
ALLOCATE(APB%rValues(SIZE(A%rValues,1),p))
ALLOCATE(APB%distance(p))
APB%distance=distance(1:p);

!HPF$ INDEPENDENT
DO i=1,p
   ! the rationale for the this 
   ! is to avoid copying padded elements


      m1=MAX(1,1-distance(i))
      m2=MIN(nrows,ncol-distance(i))

     
   SELECT CASE (dummy(i,0))
   CASE(0) ! both matrices have the diagonal
      CALL copy(A%rValues(m1:m2,dummy(i,1)),APB%rValues(m1:m2,i))
      CALL axpy(B%rValues(m1:m2,dummy(i,2)),APB%rValues(m1:m2,i))
   CASE(1) ! A has it only
      CALL copy(A%rValues(m1:m2,dummy(i,1)),APB%rValues(m1:m2,i))
   CASE(2)
      CALL copy(B%rValues(m1:m2,dummy(i,2)),APB%rValues(m1:m2,i))
   END SELECT
END DO
IF(PRESENT(freeA).AND.freeA) CALL FREEUP(A)
IF(PRESENT(freeB).AND.freeB) CALL FREEUP(B)
DEALLOCATE(dummy,distance,dist1,dist2)


END FUNCTION APB_DIA

FUNCTION APB_CRS(C,B,freeA,freeB) RESULT(APB)
  
  ! returns a pointer (APB) to the CSR representation of the sum 
  ! of two CSR sparse matrixes C and B APB=C+B 
  TYPE (CRS), POINTER :: APB
  TYPE (CRS),POINTER :: C,B
  LOGICAL, INTENT (IN), OPTIONAL :: freeA,freeB

  ! local storage
  REAL (KIND=Precision), ALLOCATABLE, DIMENSION (:) :: a
  INTEGER, ALLOCATABLE, DIMENSION (:) :: ia,ja,iw
  INTEGER :: ncol,nrow,ndummy,ierr

!0.0 check if one is null
  IF(.NOT.(ASSOCIATED(C))) THEN ! C is zero
     IF(.NOT.ASSOCIATED(B)) THEN
        APB=>NULL() ! sum is zero
        RETURN
     ELSE ! C is zero, B is not
        ALLOCATE(APB)
        APB%nNonZeros=B%nNonZeros

!!$
        ALLOCATE(APB%rValues(SIZE(B%rValues)),APB%columns(SIZE(B%columns)))
        ALLOCATE(APB%rowIndex(SIZE(B%rowIndex)))
        APB%rValues=B%rValues
        APB%columns=B%columns
        APB%rowIndex=B%RowIndex

        IF(PRESENT(freeB).AND.freeB) CALL FREEUP(B)
        RETURN
     ENDIF
  ELSE ! C is non zero
     IF(.NOT.ASSOCIATED(B)) THEN !C is non zero, B is zero
        ALLOCATE(APB)
        APB%nNonZeros=C%nNonZeros

!!$
        ALLOCATE(APB%rValues(SIZE(C%rValues)),APB%columns(SIZE(C%columns)))
        ALLOCATE(APB%rowIndex(SIZE(C%rowIndex)))
        APB%rValues=C%rValues
        APB%columns=C%columns
        APB%rowIndex=C%RowIndex
        IF(PRESENT(freeA).AND.freeA) CALL FREEUP(C)
        RETURN
     ENDIF
  ENDIF

  IF((SIZE(C%rowIndex)/=SIZE(B%rowIndex)).OR. &
      (C%nNonZeros/=B%nNonZeros)) THEN
     PRINT*,'Matrixes must have same ncol and nrow'
     PRINT*,'to be added'
     PRINT*,'FIRST MATRIX'
     PRINT*,SIZE(C%rowIndex),'X',C%nnonZeros
     PRINT*,'SECOND MATRIX'
     PRINT*,SIZE(B%rowIndex),'X',B%nnonZeros
     
     PRINT*,'Fix and restart'
  ENDIF
  ndummy=B%nNonZeros*20
  nrow=size(C%rowIndex)-1 ! rows in A

  ncol=B%nNonZeros
  ierr=1

  ALLOCATE(iw(ncol),ia(size(C%rowIndex)))
  DO WHILE(ierr/=0)

     ALLOCATE(A(ndummy),ja(ndummy))

     CALL APLB(nrow,ncol,1,C%rValues,C%columns,C%rowIndex,B%rValues,B%columns,B%rowINdex,a,ja,ia,ndummy,iw,ierr)

     IF(ierr/=0) THEN
        DEALLOCATE(A,ja)
        ndummy=ndummy+ B%nNonZeros*10! increase the size until right...
     
     ENDIF

  ENDDO
  
  ALLOCATE(APB)
  APB%nNonZeros=B%nNonZeros
!  dummy%nNonZeros=B%nNonZeros
  ndummy=IA(size(IA))-1;
!!$
  ALLOCATE(APB%rValues(ndummy),APB%columns(ndummy))
  ALLOCATE(APB%rowIndex(SIZE(IA)))

  APB%rowIndex=IA
  APB%rvalues=A(1:IA(SIZE(IA))-1);
  APB%columns=JA(1:IA(SIZE(IA))-1);

  
!!$
  

  DEALLOCATE(iw,A,ja,ia)

IF(PRESENT(freeA).AND.freeA) CALL FREEUP(C)
IF(PRESENT(freeB).AND.freeB) CALL FREEUP(B)

END FUNCTION APB_CRS

FUNCTION APB_MP(A,B,freeA,freeB)

TYPE (M_P) :: A,B,APB_MP
LOGICAL, INTENT(IN), OPTIONAL :: freeA,freeB



APB_MP%crsptr=>APB_CRS(A%crsptr,B%crsptr,freeA,freeB)
APB_MP%diaptr=>APB_DIA(A%diaptr,B%diaptr,freeA,freeB)


END FUNCTION APB_MP
  
SUBROUTINE DIAGDA(diag,A)
! performs the product diag*A in place
! wrapper to the appropriate SPARSKIT subroutine
  TYPE (CRS), INTENT (IN OUT) :: A
  REAL (KIND=Precision), DIMENSION (SIZE(A%rowIndex)-1,1,1) :: diag

  CALL DIAMUA(SIZE(A%rowIndex)-1,0,A%rValues,A%columns,A%rowIndex,diag,A%rValues,A%columns,A%rowIndex)

END SUBROUTINE DIAGDA

FUNCTION gijdj(Info,i,j,s)
TYPE (Elli) :: Info
TYPE (CRS), POINTER :: gijdj
REAL (KIND=Precision) :: dj,x,y,z,xi,eta,zeta
INTEGER, INTENT (IN) :: i,j
INTEGER :: nx,ny,nz
CHARACTER(LEN=1), INTENT (IN) :: s


! returns strip(dir=(NOR(i,j))) c2e(dir=i) gij e2c(dir=j) gc2e(dir=j) if i/=j or
! strip(other two directions) gii gc2e(dir=i) if i==j
! s controls the type of matrix it acts upon. 's'=='s' then only interior points
! are used, otherwise ghost points are assumed to be present.
! local definition

TYPE (CRS), POINTER :: A=>NULL(),B=>NULL(),C=>NULL()
REAL (KIND=Precision) :: diag
REAL (KIND=Precision), POINTER, DIMENSION (:,:,:) :: dummy
INTEGER :: idum,k,d(3,3),ii,jj
CHARACTER(LEN=5) :: type
dj=Info%dX(j)
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)


   IF(s=='b') THEN
!!$      
      IF(.NOT.g_mask(i,j)) THEN ! gij is zero, so things are easy
         ALLOCATE(dummy(nx+2,ny+2,nz+2))
         B=>operation('strip',dummy,j)
         DEALLOCATE(dummy)

         ALLOCATE(dummy(nx+2-2*d(1,j),ny+2-2*d(2,j),nz+2-2*d(3,j)))
         A=>operation('zeroi',dummy,i)
         DEALLOCATE(dummy)
         C=>ADB(A,B)
         CALL FREEUPCRS(A);CALL FREEUPCRS(B)
         k=1 ! find out the third direction
         DO WHILE(k==i.OR.k==j)
            k=k+1
         ENDDO
         ALLOCATE(dummy(nx+2-d(1,i)-2*d(1,j),ny+2-d(2,i)-2*d(2,j),nz+2-d(3,i)-2*d(3,j)))
         A=>operation('strip',dummy,k)
         gijdj=>ADB(A,C)
         DEALLOCATE(dummy)
         CALL FREEUPCRS(A);CALL FREEUPCRS(C)

         RETURN
      ENDIF
         
   
      IF(i/=j) THEN ! off diagonal elements
         ! generate gc2e along dir j
         ALLOCATE(dummy(nx+2,ny+2,nz+2))
         A=>operation('gc2eb',dummy,j,delta=dj)
         DEALLOCATE(dummy)

      ! generated e2c along same direction
         ALLOCATE(dummy(nx+2-d(1,j),ny+2-d(2,j),nz+2-d(3,j)))
         B=>operation('e2cs ',dummy,j)
         ! C=B*A
         C=>ADB(B,A)
         CALL FREEUPCRS(B);CALL FREEUPCRS(A)
         DEALLOCATE(dummy)
      ELSE ! no need for second step step
         ! generate gc2e along dir j
         ALLOCATE(dummy(nx+2,ny+2,nz+2))
         C=>operation('gc2eb',dummy,j,delta=dj)

         DEALLOCATE(dummy)



      ENDIF
      
      


! multiply by gij
! generated gij
         IF(i/=j) THEN ! needed because all metrics are defined on cell center
            ALLOCATE(dummy(d(1,j):nx+1-d(1,j),d(2,j):ny+1-d(2,j),d(3,j):nz+1-d(3,j)))
            FORALL (k=d(3,j):nz+1-d(3,j),jj=d(2,j):ny+1-d(2,j),ii=d(1,j):nx+1-d(1,j))
               dummy(ii,jj,k)=g(ii,jj,k,i,j)
            END FORALL
         ELSE
            ALLOCATE(dummy(d(1,j):nx+1,d(2,j):ny+1,d(3,j):nz+1))
            FORALL(k=d(3,j):nz+1,jj=d(2,j):ny+1,ii=d(1,j):nx+1)
               dummy(ii,jj,k)=g(ii,jj,k,j,j)
            END FORALL
                     
         ENDIF


         CALL DIAGDA(dummy,C)

         
   ! IF i==j we need to strip other two dimensions. 
   ! IF i/=j then we need to c2e along(i) and strip along remaining direction

   ! strip along k/=(i,j). First find out which direction it is
         k=1
         DO WHILE(k==i.OR.k==j)
            k=k+1
         ENDDO

         
         ! dummy has the right dimension
         B=>operation('strip',dummy,k)

         A=>ADB(B,C)

         CALL FREEUPCRS(C);

         CALL FREEUPCRS(B)

         DEALLOCATE(dummy)

         IF(i==j) THEN ! apply the last strip
            ALLOCATE(dummy(d(1,j)+d(1,k):nx+1-d(1,k), &
                           d(2,j)+d(2,k):ny+1-d(2,k), &
                           d(3,j)+d(3,k):nz+1-d(3,k)))
            k=k+1
            DO WHILE(k==i)
               k=k+1
            END DO

            B=>operation('strip',dummy,k)

            gijdj=>ADB(B,A)
            DEALLOCATE(dummy)
         ELSE

            ALLOCATE(dummy(d(1,j)+d(1,k):nx+1-d(1,j)-d(1,k), &
                           d(2,j)+d(2,k):ny+1-d(2,j)-d(2,k), &
                           d(3,j)+d(3,k):nz+1-d(3,j)-d(3,k)))
            
            
            B=>operation('c2eb ',dummy,i)

            gijdj=>ADB(B,A)
            DEALLOCATE(dummy)
         ENDIF

         ! release everthing else

         

      ELSE
         PRINT*,'option s not yet implemented in gijdj'
         PRINT*,'use G_(Info) instead to construct the'
         PRINT*,'gradient operator'
         STOP
      ENDIF

      CALL FREEUPCRS(A);CALL FREEUPCRS(B);CALL FREEUPCRS(C)

    END FUNCTION gijdj
   
    SUBROUTINE FREEUPCRS(A)
! releases memory taken by a CRS pointer
      TYPE (CRS),POINTER :: A
      IF(ASSOCIATED(A)) THEN
         IF(ASSOCIATED(A%rValues)) THEN
            DEALLOCATE(A%rValues); nullify (A%rValues)
         ENDIF
         IF(ASSOCIATED(A%columns)) THEN
            DEALLOCATE(A%columns); nullify (A%columns)
         ENDIF
         IF(ASSOCIATED(A%rowIndex)) THEN
            DEALLOCATE(A%rowIndex); nullify (A%rowIndex)
         ENDIF
         DEALLOCATE(A);NULLIFY(A) ! this releases the memory taken by
                                  ! the metadata and statically allocated
                                  ! entries of the structure
      ENDIF
    END SUBROUTINE FREEUPCRS

    SUBROUTINE FREEUPDIA(A)
! releases memory taken by a DIA pointer
      TYPE (DIA),POINTER :: A
      IF(ASSOCIATED(A)) THEN
         IF(ASSOCIATED(A%rValues)) THEN
            DEALLOCATE(A%rValues); nullify (A%rValues)
         ENDIF
         IF(ASSOCIATED(A%distance)) THEN
            DEALLOCATE(A%distance); nullify (A%distance)
         ENDIF
         DEALLOCATE(A);NULLIFY(A) ! this releases the memory taken by
                                  ! the metadata and statically allocated
                                  ! entries of the structure
      ENDIF
    END SUBROUTINE FREEUPDIA
    SUBROUTINE FREEUPM_P(A)
      ! releases memory taken by a M_P pointer; note that the input is not a pointer
      TYPE(M_P) :: A

         IF(ASSOCIATED(A%crsptr)) THEN
            CALL FREEUPCRS(A%crsptr)
         ENDIF
         IF(ASSOCIATED(A%diaptr)) THEN
            CALL FREEUPDIA(A%diaptr)
         ENDIF
!!$         DEALLOCATE(A)

    END SUBROUTINE FREEUPM_P

SUBROUTINE CRS2DIA(G,maxdia,stats)

TYPE (M_P) :: G
INTEGER, OPTIONAL, INTENT(IN) :: maxdia
LOGICAL, OPTIONAL, INTENT(IN) :: stats
TYPE (CRS), POINTER :: crsptr
TYPE (DIA), POINTER :: diaptr
INTEGER :: n,job,idiag,ndiag
REAL (KIND=Precision), DIMENSION (:), POINTER :: rValues
REAL (KIND=Precision), DIMENSION (:,:), POINTER :: diag
INTEGER, DIMENSION (:), POINTER :: ind,distance,columns,rowIndex
LOGICAL :: PRINT_STATS

! from a matrix stored entirely as a CRS, extracts the elements can can be
! efficiently stored in DIA format, leaving behind the possible scattered 
! one in the CRS portion

INTERFACE 
   subroutine csrdia (n,idiag,job,a,ja,ia,ndiag, &
                        diag,distance,ao,jao,iao,ind)
      real*8 diag(ndiag,idiag), a(*), ao(*)
      integer ia(*), ind(*), ja(*), jao(*), iao(*), distance(*)
    end subroutine csrdia
 END INTERFACE

IF(.NOT.PRESENT(stats)) THEN
   PRINT_STATS=.FALSE.
ELSE
   PRINT_STATS=.TRUE.
ENDIF

! assign local pointers
IF(.NOT.ASSOCIATED(G%diaptr)) THEN
   ALLOCATE(G%diaptr)
ELSE ! discard previous content
   CALL FREEUP(G%diaptr)
   ALLOCATE(G%diaptr)
ENDIF
crsptr=>G%crsptr
diaptr=>G%diaptr
IF(PRINT_STATS) THEN
   PRINT*,'prior to rebalancing'
   PRINT*,'the crs is as such'
   PRINT*,'size rvalues'
   PRINT*,SIZE(crsptr%rValues)
   PRINT*,'size columns'
   PRINT*,SIZE(crsptr%columns)
ENDIF

n=SIZE(crsptr%rowIndex)-1;
ALLOCATE(ind(2*n-1))
IF(.NOT.PRESENT(maxdia)) THEN
   idiag=20 ! this should be enough for the worst possible case
ELSE
   idiag=maxdia
ENDIF
job=11 ! select diagonals automatically and leaves the rest in the CRS part
ALLOCATE(distance(idiag))
ALLOCATE(diag(n,idiag))

CALL csrdia(n,idiag,job,crsptr%rvalues,crsptr%columns,crsptr%RowIndex,&
     n,diag,distance,crsptr%rvalues,crsptr%columns,crsptr%rowIndex,ind)
DEALLOCATE(ind)
! save results in DIA 

ALLOCATE(diaptr%rValues(n,idiag))
diaptr%rValues=diag(1:n,1:idiag)
DEALLOCATE(diag)

ALLOCATE(diaptr%distance(idiag))
diaptr%distance=distance(1:idiag)
DEALLOCATE(distance)
diaptr%ncol=G%crsptr%nNonZeros
! resize CRS part

NULLIFY(crsptr)

G%crsptr=>REALLOCATE(G%crsptr)

crsptr=>G%crsptr
IF(PRINT_STATS) THEN
   IF(ASSOCIATED(crsptr)) THEN
      PRINT*,'rebalancing completed'
      PRINT*,'now the crs is as such'
      PRINT*,'size rvalues'
      PRINT*,SIZE(crsptr%rValues)
      PRINT*,'size columns'
      PRINT*,SIZE(crsptr%columns)
      PRINT*,'size diags'
      PRINT*,SIZE(diaptr%distance)
   ELSE
      PRINT*,'the crs is empty'
   ENDIF
      PRINT*,'diagonals extracted'
      PRINT*,diaptr%distance

ENDIF

END SUBROUTINE CRS2DIA

FUNCTION REALLOCATE_CRS(G)

TYPE(CRS), POINTER :: G,REALLOCATE_CRS

INTEGER :: n
! used to shrink a crs structure. Note that G is
! deallocated on exit

n=G%rowIndex(SIZE(G%rowIndex))-1
IF(n>0) THEN 
   ALLOCATE(REALLOCATE_CRS)

   ALLOCATE(REALLOCATE_CRS%rValues(n))
   ALLOCATE(REALLOCATE_CRS%columns(n))
   ALLOCATE(REALLOCATE_CRS%rowIndex(SIZE(G%rowIndex)))
   
   REALLOCATE_CRS%rowIndex=G%rowIndex
   
   REALLOCATE_CRS%columns=G%columns(1:n)
   REALLOCATE_CRS%rValues=G%rValues(1:n)
   
   REALLOCATE_CRS%nNonZeros=G%nNOnzeros
ELSE ! for an empty matrix returns a null pointer
   REALLOCATE_CRS=>NULL()
ENDIF

CALL FREEUP(G)

END FUNCTION REALLOCATE_CRS



    FUNCTION EXTRAPOLATE_BC(Info)
! returns the pointer to the CRS representation of the operator that
! acts on p on the interior points and return p on the interior+ghost points.
! The values on the ghost points are establihed as follows
! on boundary i if Info%bc(i)%bc_type==DR then ghostpoint=0 (since the extrapolation is given)
! on boundary i if Info%bc(i)%bc_type==NM then ghostpoint=depends on interior layers
! the matrix is (nx+2)*(ny+2)*(nz+2) X nx*ny*nz
! interface declaration
TYPE (Elli) :: Info
TYPE (CRS), POINTER :: EXTRAPOLATE_BC
! local declarations
TYPE (CRS), POINTER :: A=>NULL(),B=>NULL(),C=>NULL()
REAL (KIND=Precision), DIMENSION (:,:,:), POINTER :: dummy
INTEGER :: nx,ny,nz
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
ALLOCATE(dummy(nx,ny,nz))

IF(Info%nDim==1) THEN
   EXTRAPOLATE_BC=>BoC(Info,1)
ELSE
   A=>BoC(Info,1)
   
   B=>BoC(Info,2)
   
   IF(Info%nDim==2) THEN
      
      EXTRAPOLATE_BC=>APB(A,B)
      CALL FREEUPCRS(A);CALL FREEUPCRS(B)
   ELSE
      C=>APB(A,B)
      
      CALL FREEUPCRS(A)
      CALL FREEUPCRS(B)      

   

   
   

      A=>BoC(Info,3)
      
      EXTRAPOLATE_BC=>APB(A,C)
      CALL FREEUPCRS(C)
      CALL FREEUPCRS(A)

      
      
   ENDIF
ENDIF
DEALLOCATE(dummy)
END FUNCTION EXTRAPOLATE_BC

FUNCTION BoC(Info,d)
! CRS representation of the matrix extrapolating the pressure from internal points to internal+ghost
! along direction i
! interface declaration
TYPE (Elli) :: Info
INTEGER :: d
TYPE (CRS), POINTER :: BoC
! local declaration
INTEGER nx,ny,nz,boundary(2),N,nNonZeros,i,j,k,l,m
REAL (KIND=Precision) :: dx,dy,dz,fl(4),betax,betay,betaz
REAL (KIND=Precision),  DIMENSION (:), POINTER :: rValues
INTEGER, DIMENSION (:), POINTER :: columns

dx=Info%dX(1);dy=Info%dX(2);dz=Info%dX(3)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
fl=zero
boundary=Info%bc((d-1)*2+1:d*2)%bc_type

! determine size of matrix and allocate space

ALLOCATE(BoC)

BoC%nNonZeros=(nx)*(ny)*(nz) ! ncol of matrix
N=0 ! 
! size of matrix depends on three cases 
IF(ALL(boundary==DR)) THEN
   nNonZeros=N
ELSEIF(ALL(boundary==NM)) THEN
   SELECT CASE (d)
   CASE (1)
      nNonZeros=N+10*(ny+2)*(nz+2)
   CASE (2)
      nNonZeros=N+10*(nx+2)*(nz+2)
   CASE (3)
      nNonZeros=N+10*(ny+2)*(nx+2)
   CASE DEFAULT
      PRINT*,'BC called with', i,' out of bounds [1,3]'
      STOP

   END SELECT
   
ELSE ! one DR and one NM
   SELECT CASE (d)
   CASE (1)
      nNonZeros=N+5*(ny+2)*(nz+2)
   CASE (2)
      nNonZeros=N+5*(nx+2)*(nz+2)
   CASE (3)
      nNonZeros=N+5*(ny+2)*(nx+2)
   CASE DEFAULT
      PRINT*,'BC called with', i,' out of bounds [1,3]'
      STOP

   END SELECT
   
ENDIF

N=(nx+2)*(ny+2)*(nz+2) ! nrow
ALLOCATE(BoC%columns(nNonZeros),BoC%rValues(nNonZeros))
ALLOCATE(BoC%rowIndex(N+1))
columns=>BoC%columns;rValues=>BoC%rValues

nNonZeros=1
! fill in the matrix. There are 12 combinations

   IF(ALL(boundary==DR)) THEN ! there is no need to do anything. the output is a zero matrix
      
      rvalues=one;columns=1;BoC%rowIndex=nNonZeros ! note, rvalues and columns are an empty matrix
      RETURN
   ELSE ! need to distinguish



      SELECT CASE (d)
         
      CASE (1) ! boundaries in x direction
         DO i=1,2*(Info%ndim-1)
            IF(Info%bc(i+2)%bc_type==NM) THEN
               fl(i)=one ! need to add p(:,:,1) to the extrapolation
            ELSE
               fl(i)=zero
            ENDIF
         ENDDO

         IF(ALL(boundary-(/NM,DR/)==0)) THEN



            DO k=0,nz+1
               DO j=0,ny+1
                  i=0
                  l=(k-1)*nx*ny+(j-1)*nx+i+1 ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  
                  betax=0.25d0*(g(i,j,k,1,2)+g(i+1,j,k,1,2))/g(i+1,j,k,1,1)/dy*dx
                  IF(Info%nDim>2) THEN
                     betay=0.25d0*(g(i,j,k,1,3)+g(i+1,j,k,1,3))/g(i+1,j,k,1,1)/dz*dx
                     CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,j,k,ny,nz,nx,nx*ny)
                  ELSEIF(Info%nDim>1) THEN
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,j,ny,nx)
                  ELSE
                     IF(k==1.AND.j==1) THEN
                        rValues(nNonZeros:nNonzeros)=one
                        columns(nNonzeros:nNonzeros)=(/l/)
                        nNonzeros=nNonzeros+1
                     ENDIF
                  ENDIF




               
               DO i=1,nx+1
                     l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                     m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                     BoC%rowIndex(m)=nNonZeros

                  ENDDO
               ENDDO
            ENDDO
            
            
            
         ELSEIF(ALL(boundary-(/DR,NM/)==0)) THEN
            DO k=0,nz+1
               DO j=0,ny+1
                  DO i=0,nx
                     l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                     m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                     BoC%rowIndex(m)=nNonZeros
                  ENDDO

               
                  i=nx+1
                  l=(k-1)*nx*ny+(j-1)*nx+i-1 ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  
                  betax=-0.25d0*(g(i-1,j,k,1,2)+g(i,j,k,1,2))/g(i,j,k,1,1)/dy*dx
                  IF(Info%nDim>2) THEN
                     betay=-0.25d0*(g(i-1,j,k,1,3)+g(i,j,k,1,3))/g(i,j,k,1,1)/dz*dx
                     CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,j,k,ny,nz,nx,nx*ny)
                  ELSEIF(Info%nDim>1) THEN
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,j,ny,nx)
                     betay=zero
                  ELSE
                     IF(k==1.AND.j==1) THEN
                        rValues(nNonZeros:nNonzeros)=one
                        columns(nNonzeros:nNonzeros)=(/l/)
                        nNonzeros=nNonzeros+1
                     ENDIF
                        
                  ENDIF


               ENDDO
            ENDDO

         ELSE 

            DO k=0,nz+1
               DO j=0,ny+1
                  i=0
                  l=(k-1)*nx*ny+(j-1)*nx+i+1 ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  
                  betax=0.25d0*(g(i,j,k,1,2)+g(i+1,j,k,1,2))/g(i+1,j,k,1,1)/dy*dx
                  IF(Info%nDim>2) THEN
                     betay=0.25d0*(g(i,j,k,1,3)+g(i+1,j,k,1,3))/g(i+1,j,k,1,1)/dz*dx
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,j,k,ny,nz,nx,nx*ny)
                  ELSEIF(Info%nDim>1) THEN
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,j,ny,nx)
                     betay=zero
                  ELSE
                     IF(k==1.AND.j==1) THEN
                        rValues(nNonZeros:nNonzeros)=one
                        columns(nNonzeros:nNonzeros)=(/l/)
                        nNonzeros=nNonzeros+1
                     ENDIF

                  ENDIF


                  
                  
                  DO i=1,nx
                     l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                     m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                     BoC%rowIndex(m)=nNonZeros
                  ENDDO
                  i=nx+1
                  l=(k-1)*nx*ny+(j-1)*nx+i-1 ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  
                  betax=-0.25d0*(g(i-1,j,k,1,2)+g(i,j,k,1,2))/g(i,j,k,1,1)/dy*dx
                  IF(Info%nDim>2) THEN
                     betay=-0.25d0*(g(i-1,j,k,1,3)+g(i,j,k,1,3))/g(i,j,k,1,1)/dz*dx
                     CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,j,k,ny,nz,nx,nx*ny)
                  ELSEIF(Info%nDim>1) THEN
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,j,ny,nx)
                     betay=zero
                  ELSE
                     IF(k==1.AND.j==1) THEN
                        rValues(nNonZeros:nNonzeros)=one
                        columns(nNonzeros:nNonzeros)=(/l/)
                        nNonzeros=nNonzeros+1
                     ENDIF

                  ENDIF


               ENDDO
            ENDDO
            

         END IF
         
      CASE (2) ! boundaries in y direction
         DO i=1,2
            IF(Info%bc(i)%bc_type==NM) THEN
               fl(i)=one ! need to add p(:,:,1) to the extrapolation
            ELSE
               fl(i)=zero
            ENDIF
         ENDDO
         IF(Info%nDim>2) THEN
            DO i=3,4
               IF(Info%bc(i+2)%bc_type==NM) THEN
                  fl(i)=one ! need to add p(:,:,1) to the extrapolation
               ELSE
                  fl(i)=zero
               ENDIF
            ENDDO
         ENDIF

         IF(ALL(boundary-(/NM,DR/)==0)) THEN

            DO k=0,nz+1
               j=0
               DO i=0,nx+1
                  l=(k-1)*nx*ny+(j)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          

                  betax=0.25d0*(g(i,j,k,2,1)+g(i,j+1,k,2,1))/g(i,j+1,k,2,2)/dx*dy
                  IF(Info%nDim>2) THEN
                     betay=0.25d0*(g(i,j,k,2,3)+g(i,j+1,k,2,3))/g(i,j+1,k,2,2)/dz*dy
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,k,nx,nz,1,nx*ny)
                  ELSE
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,i,nx,1)
                     betay=zero
                  ENDIF

               ENDDO
               DO j=1,ny+1
                  DO i=0,nx+1
                  l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                     
                  ENDDO
               ENDDO
            ENDDO
            
            
         ELSEIF(ALL(boundary-(/DR,NM/)==0)) THEN
            DO k=0,nz+1
               DO j=0,ny
                  DO i=0,nx+1
                  l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                     
                  ENDDO
               ENDDO
               j=ny+1
               DO i=0,nx+1
                  l=(k-1)*nx*ny+(j-1-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  betax=-0.25d0*(g(i,j,k,2,1)+g(i,j-1,k,2,1))/g(i,j,k,2,2)/dx*dy
                  IF(Info%nDim>2) THEN
                     betay=-0.25d0*(g(i,j,k,2,3)+g(i,j-1,k,2,3))/g(i,j,k,2,2)/dz*dy
                     CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,k,nx,nz,1,nx*ny)
                  ELSE
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,i,nx,1)
                     betay=zero
                  ENDIF

                  
                  
               ENDDO
            ENDDO
            
         ELSE
            DO k=0,nz+1
               j=0
               DO i=0,nx+1
                  l=(k-1)*nx*ny+(j+1-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  betax=0.25d0*(g(i,j,k,2,1)+g(i,j+1,k,2,1))/g(i,j+1,k,2,2)/dx*dy
                  IF(Info%nDim>2) THEN
                     betay=0.25d0*(g(i,j,k,2,3)+g(i,j+1,k,2,3))/g(i,j+1,k,2,2)/dz*dy
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,k,nx,nz,1,nx*ny)
                  ELSE
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,i,nx,1)
                     betay=zero
                  ENDIF

                  

               ENDDO
               
               DO j=1,ny
                  DO i=0,nx+1
                  l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                     
                  ENDDO
               ENDDO
               
               j=ny+1
               DO i=0,nx+1
                  l=(k-1)*nx*ny+(j-1-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  betax=-0.25d0*(g(i,j,k,2,1)+g(i,j-1,k,2,1))/g(i,j,k,2,2)/dx*dy
                  IF(Info%nDim>2) THEN
                     betay=-0.25d0*(g(i,j,k,2,3)+g(i,j-1,k,2,3))/g(i,j,k,2,2)/dz*dy
                     CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,k,nx,nz,1,nx*ny)
                  ELSE
                     IF(k==1) &
                          CALL LOOP2D(rValues,columns,betax,fl,nNonZeros,l,i,nx,1)
                     betay=zero
                  ENDIF

                  
                  
               ENDDO
            ENDDO
            
         ENDIF
      CASE (3) ! boundaries in z direction
         
         DO i=1,4
            IF(Info%bc(i)%bc_type==NM) THEN
               fl(i)=one ! need to add p(:,:,1) to the extrapolation
            ELSE
               fl(i)=zero
            ENDIF
         ENDDO

         IF(ALL(boundary-(/NM,DR/)==0)) THEN



            k=0
            DO j=0,ny+1
               DO i=0,nx+1
                  l=(k+1-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          
                  betax=0.25d0*(g(i,j,k,3,1)+g(i,j,k+1,3,1))/g(i,j,k+1,3,3)/dx*dz
                  betay=0.25d0*(g(i,j,k,3,2)+g(i,j,k+1,3,2))/g(i,j,k+1,3,3)/dy*dz
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,j,nx,ny,1,nx)
               ENDDO
            ENDDO
            DO k=1,nz+1
               DO j=0,ny+1
                  DO i=0,nx+1
                     l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                     m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                     BoC%rowIndex(m)=nNonZeros
                  ENDDO
               ENDDO
            ENDDO

                     

         ELSEIF(ALL(boundary-(/DR,NM/)==0)) THEN
            DO k=0,nz
               DO j=0,ny+1
                  DO i=0,nx+1
                     l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                     m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                     BoC%rowIndex(m)=nNonZeros
                  ENDDO
               ENDDO
            ENDDO

            k=nz+1
            DO j=0,ny+1
               DO i=0,nx+1
                  l=(k-1-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          

                  betax=-0.25d0*(g(i,j,k-1,3,1)+g(i,j,k,3,1))/g(i,j,k,3,3)/dx*dz
                  betay=-0.25d0*(g(i,j,k-1,3,2)+g(i,j,k,3,2))/g(i,j,k,3,3)/dy*dz
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,j,nx,ny,1,nx)
                  
               ENDDO
            ENDDO

            
         ELSE 

            k=0
            DO j=0,ny+1
               DO i=0,nx+1
                  l=(k+1-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          

                  betax=0.25d0*(g(i,j,k,3,1)+g(i,j,k+1,3,1))/g(i,j,k+1,3,3)/dx*dz
                  betay=0.25d0*(g(i,j,k,3,2)+g(i,j,k+1,3,2))/g(i,j,k+1,3,3)/dy*dz
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,j,nx,ny,1,nx)
                  
               ENDDO
            ENDDO


            DO k=1,nz
               DO j=0,ny+1
                  DO i=0,nx+1
                     l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                     m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                     BoC%rowIndex(m)=nNonZeros

                  ENDDO
               ENDDO
            ENDDO

            k=nz+1
            DO j=0,ny+1
               DO i=0,nx+1
                  l=(k-1-1)*nx*ny+(j-1)*nx+i ! vectorized column index
                  m=k*(nx+2)*(ny+2)+j*(nx+2)+i+1 ! vectorized row index
                  BoC%rowIndex(m)=nNonzeros          

                  betax=-0.25d0*(g(i,j,k-1,3,1)+g(i,j,k,3,1))/g(i,j,k,3,3)/dx*dz
                  betay=-0.25d0*(g(i,j,k-1,3,2)+g(i,j,k,3,2))/g(i,j,k,3,3)/dy*dz
                  CALL LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,j,nx,ny,1,nx)
                  
               ENDDO
            ENDDO


         END IF

      END SELECT
   END IF
   BoC%rowIndex(m+1)=nNonZeros
   BoC=>FILTER_Z_P(BoC,1.d-14) ! useful in 2D 
 END FUNCTION BoC

SUBROUTINE LOOP(rValues,columns,betax,betay,fl,nNonZeros,l,i,j,nx,ny,stx,sty)
! fills in the nonzeros values in a row 
REAL (KIND=Precision), DIMENSION (:), INTENT (IN OUT) :: rvalues
INTEGER, DIMENSION (:), INTENT (IN OUT) :: columns
INTEGER, INTENT (IN) :: l,i,j,nx,ny,stx,sty
INTEGER, INTENT (IN OUT) :: nNonZeros
REAL (KIND=Precision), INTENT (IN) :: betax,betay,fl(4)
INTEGER :: nlx,nly
nlx=max(nx,2);nly=max(ny,2) ! hack to avoid trouble if nx or ny ==1
                            
IF(MOD(j,ny+1)*MOD(i,nx+1)/=0) THEN
   IF(MOD(i-1,nlx-1)+MOD(j-1,nly-1)==0) THEN ! vertices
      IF(i==1.AND.j==1) THEN
         rValues(nNonZeros:nNonZeros+2)=(/one-fl(1)*betax-fl(3)*betay,betax,betay/)
         columns(nNonZeros:nNonZeros+2)=(/l,l+stx,l+sty/)
         nNonZeros=nNonZeros+3
      ELSEIF(i==1.AND.j==ny) THEN
         rValues(nNonZeros:nNonZeros+2)=(/-betay,one-fl(1)*betax+fl(4)*betay,betax/)
         columns(nNonZeros:nNonZeros+2)=(/l-sty,l,l+stx/)
         nNonZeros=nNonZeros+3
      ELSEIF(i==nx.AND.j==1) THEN
         rValues(nNonZeros:nNonZeros+2)=(/-betax,one+fl(2)*betax-fl(3)*betay,betay/)
         columns(nNonZeros:nNonZeros+2)=(/l-stx,l,l+sty/)
         nNonZeros=nNonZeros+3
      ELSE
         rValues(nNonZeros:nNonZeros+2)=(/-betay,-betax,one+fl(2)*betax+fl(4)*betay/)
         columns(nNonZeros:nNonZeros+2)=(/l-sty,l-stx,l/)
         nNonZeros=nNonZeros+3
      ENDIF
   ELSEIF(mod(i-1,nlx-1)==0) THEN
      IF(i==1) THEN
         rValues(nNonZeros:nNonZeros+3)=(/-betay,one-fl(1)*betax,betax,betay/)
         columns(nNonZeros:nNonZeros+3)=(/l-sty,l,l+stx,l+sty/)
         nNonZeros=nNonZeros+4
      ELSE
         rValues(nNonZeros:nNonZeros+3)=(/-betay,-betax,one+fl(2)*betax,betay/)
         columns(nNonZeros:nNonZeros+3)=(/l-sty,l-stx,l,l+sty/)
         nNonZeros=nNonZeros+4
      ENDIF
   ELSEIF(mod(j-1,nly-1)==0) THEN
      IF(j==1) THEN
         rValues(nNonZeros:nNonZeros+3)=(/-betax,one-fl(3)*betay,betax,betay/)
         columns(nNonZeros:nNonZeros+3)=(/l-stx,l,l+stx,l+sty/)
         nNonZeros=nNonZeros+4
      ELSE
         rValues(nNonZeros:nNonZeros+3)=(/-betay,-betax,one+fl(4)*betay,betax/)
         columns(nNonZeros:nNonZeros+3)=(/l-sty,l-stx,l,l+stx/)
         nNonZeros=nNonZeros+4
      ENDIF
   ELSE
      rValues(nNonZeros:nNonZeros+4)=(/-betay,-betax,one,betax,betay/)
      columns(nNonZeros:nNonZeros+4)=(/l-sty,l-stx,l,l+stx,l+sty/)
      nNonZeros=nNonZeros+5
   END IF
END IF

END SUBROUTINE LOOP

SUBROUTINE LOOP2D(rValues,columns,betax,fl,nNonZeros,l,i,nx,stx)
! fills in the nonzeros values in a row 
REAL (KIND=Precision), DIMENSION (:), INTENT (IN OUT) :: rvalues
INTEGER, DIMENSION (:), INTENT (IN OUT) :: columns
INTEGER, INTENT (IN) :: l,i,nx,stx
INTEGER, INTENT (IN OUT) :: nNonZeros
REAL (KIND=Precision), INTENT (IN) :: betax,fl(4)
INTEGER :: nlx,nly

                            
IF(MOD(i,nx+1)/=0) THEN
   IF(MOD(i-1,nx-1)==0) THEN ! vertices
      IF(i==1) THEN
         rValues(nNonZeros:nNonZeros+1)=(/one-fl(1)*betax,betax/)
         columns(nNonZeros:nNonZeros+1)=(/l,l+stx/)
         nNonZeros=nNonZeros+2
      ELSEIF(i==nx) THEN
         rValues(nNonZeros:nNonZeros+1)=(/-betax,one+fl(2)*betax/)
         columns(nNonZeros:nNonZeros+1)=(/l-stx,l/)
         nNonZeros=nNonZeros+2
      ENDIF
   ELSE
      rValues(nNonZeros:nNonZeros+2)=(/-betax,one,betax/)
      columns(nNonZeros:nNonZeros+2)=(/l-stx,l,l+stx/)
      nNonZeros=nNonZeros+3
   END IF
END IF

END SUBROUTINE LOOP2D

SUBROUTINE copy(X,Y)
! wrapper for copy
  REAL (KIND=Precision), DIMENSION (:):: X,Y
  call dcopy(SIZE(X),X,1,Y,1)

END SUBROUTINE copy

SUBROUTINE axpy(X,Y)
 ! wrapper
  REAL (KIND=Precision), DIMENSION (:):: X,Y
  CALL daxpy(SIZE(X),one,X,1,Y,1)
END SUBROUTINE axpy

!********************
END MODULE Matrix_Oper
!********************

